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Abstract. 

For  quasistatic  stress  problems  two  alternative  constitutive  relationships  expressing  the  stress  in  a 
linear  viscoelastic  solid  body  as  a  linear  functional  of  the  strain  are  derived.  In  conjunction  with 
the  equations  of  equilibrium,  these  form  the  mathematical  models  for  the  stress  problems.  These 
models  are  first  discretized  in  the  space  domain  using  a  finite  element  method  and  semi-discrete 
error  estimates  are  presented  corresponding  to  each  constitutive  relationship.  Through  the  use 
respectively  of  quadrature  rules  and  finite  difference  replacements  each  semi-discrete  scheme  is 
fully  discretized  into  the  time  domain  so  that  two  practical  algorithms  suitable  for  the  numerical 
stress  analysis  of  linear  viscoelastic  solids  are  produced.  The  semi-discrete  estimates  are  then  also 
extended  into  the  time  domain  to  give  spatially  H 1  error  estimates  for  each  algorithm. 

The  numerical  schemes  are  predicated  on  exact  analytical  solutions  for  a  simple  model  problem, 
and  finally  on  design  data  for  a  real  polymeric  material. 

1  Introduction. 

1.1  Outline. 

This  paper  is  concerned  with  the  computational  modelling  of  problems  of  solid  mechanics  in  which 
the  material  is  assumed  to  be  linear  viscoelastic  and  the  deformation  is  quasistatic  (i.e.  the  inertia 
term  in  the  momentum  equations  is  neglected,  providing  equilibrium  equations).  For  any  problem 
of  this  type  our  purpose  here  is  to  produce  a  mathematical  model,  to  discretize  this  to  produce  a 
numerical  algorithm,  to  derive  a  priori  error  estimates  for  the  error  associated  with  the  numerical 
solution  and  to  apply  the  scheme  to  a  number  of  test  problems  for  numerical  verification.  We 
have  additionally  applied  the  algorithm  to  a  problem  involving  experimental  data  from  a  specific 
polymeric  material. 

Two  forms  of  the  constitutive  relation  are  proposed,  giving  two  models.  These  involve  the 
rate  of  change  with  time  respectively  of  the  stress  relaxation  function  [1,3,4]  and  the  displacement 
vector.  Both  models  are  first  discretized  in  the  space  variables  using  a  Galerkin  technique,  thus 
setting  up  semi-discrete  approximating  problems.  These  are  further  discretized  in  time  through 
the  application  respectively  of  quadrature  rules  and  finite  difference  replacements,  thus  producing 
two  fully  discrete  approximating  formulations.  The  errors  in  these  are  analysed  and,  in  §§3,4,  a 
priori  estimates,  for  the  semi-discrete  and  fully  discrete  schemes  respectively,  are  derived.  The 
results  of  numerical  experiments  are  given  in  §5  for  some  problems  where  the  exact  solutions  have 
been  obtained  through  application  of  the  correspondence  principle  of  Schapery  [3,4],  and  for  the 
specific  polymeric  material  as  above. 

1.2  Viscoelastic  materials. 

The  stress  analysis  of  solid  materials  is  usually  restricted  to  the  cases  either  of  linear  elasticity, 
under  the  assumptions  that  the  deflections  considered  are  in  some  sense  small  and  that  the  stresses 
occuring  in  the  body  never  exceed  the  yield  stress  of  the  material,  or  of  plasticity  where  the  yield 
stress  is  exceeded  and  some  permanent  deformation  of  the  material  takes  place.  Whereas  linear 
elasticity  utilises  a  simple  relationship  (Hooke’s  law)  between  stress  and  strain,  plasticity  is  a 
nonlinear  phenomena  and  in  flow  theory  the  constitutive  relationship  is  strain  rate  dependent  and 
thus  depends  on  the  loading  history. 
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There  is,  however,  a  class  of  materials,  solids  and  fluids,  the  deformation  of  which  can  be 
modelled  using  linear  theory,  but  for  which  linear  elasticity  is  not  appropriate.  These  materials 
are  termed  viscoelastic  since  in  deformation  they  display  both  elastic  and  viscous  flow  properties. 

Perhaps  the  most  common  examples  of  viscoelastic  materials  in  engineering  science  are  the 
thermoplastic  polymers.  In  simple  terms  such  materials  may  be  considered  to  be  composed  of  many 
so  called  long  chain  molecules  arranged  more  or  less  at  random  and  intertwining  with  each  other 
in  a  “spaghetti”  like  manner.  Neighbouring  chains  may  be  joined  to  each  other  by  cross-ltnking . 
When  a  mechanical  load  is  applied  to  such  a  material  there  will  be  an  elastic  deflection ,  due  to  the 
the  ability  of  each  chain  to  stretch,  and  a  viscous  flow,  caused  by  the  sliding  of  the  molecules  over 
one  another;  the  extent  to  which  molecules  may  slide  in  this  way  is  determined  by  the  amount  of 
cross-linking.  This  dual  (elastic  plus  viscous)  effect  gives  rise  to  the  term  viscoelasticity .  If  the 
degree  of  cross-linking  is  very  small  the  flow  may  continue  unchecked  and  the  material  is  considered 
to  be  a  viscoelastic  liquid.  Conversely,  if  the  degree  is  high  the  flow  will  eventually  cease  (unless 
the  applied  load  is  large  enough  to  rupture  the  material)  and  the  material  is  a  viscoelastic  solid. 
In  practice  this  distinction  may  be  somewhat  blurred,  particularly  in  nonisothermal  contexts  near 
the  melt  temperatures  of  viscoelastic  materials. 

In  the  present  paper  we  consider  only  compressible,  isotropic  viscoelastic  media  undergoing 
“small”  deformations  subject  to  isothermal,  quasistatic  conditions. 

1.3  Linear  viscoelasticity. 

Experimental  observations  of  viscoelastic  materials  under  loads,  see  e.g.  [3],  suggest  that  at  any 
given  time  the  stress  at  a  point  in  a  viscoelastic  medium  depends  not  only  on  the  current  strain 
at  that  point  but  upon  the  entire  strain  history  from  the  instant  that  the  material  was  in  its 
original  unstressed  state.  For  this  reason  viscoelastic  materials  are  described  as  having  memory. 
Mathematically  this  suggests  that  the  stress  may  be  written  as  a  functional  of  the  strain,  or  vice- 
versa-,  a  linear  functional  gives  rise  to  the  linear  viscoelastic  model.  When  the  material  function 
involved  in  this  functional  depends  only  upon  the  relative  time  (i.e.  the  algebraic  difference  between 
the  present  time  and  the  time  in  the  deformation  or  loading  history)  the  material  so  described  is 
termed  “non-ageing”. 

In  order  to  arrive  at  an  explicit  representation  of  this  functional  we  have  the  following: 
Hypothesis  1.3.1.  The  Boltzmann  superposition  principle  [4], 

For  a  linear,  non-ageing  viscoelastic  body  the  stress  arising  at  position  x  and  time  t  from  a  sequence 
of  strain  increments  {A^x,**)}^1,  </v-i  <  t,  applied  to  the  body  may  be  expressed  as  the  sum 
of  the  stresses  {A(r‘(x,  f)}^1  induced  by  each  individual  strain  increment. 

Two  readily  observable  phenomena  characterizing  a  viscoelastic  material  are  those  of  creep 
and  stress  relaxation.  Consider  a  slim  viscoelastic  cylindrical  bar  subjected  to  a  constant  axial 
(tensile)  load.  In  the  deformation  one  would  observe  an  instantaneous  elastic  deflection  followed 
by  a  viscous  flow,  possibly  up  to  some  equilibrium  deformation;  this  is  termed  creep.  Now,  consider 
a  similar  bar  subjected  to  a  constant  axial  strain.  The  stress  in  the  bar  will  instantaneously  increase 
as  an  elastic  response  and  will  thereafter  decay  to  some  constant  (possibly  zero)  value;  this  is  stress 
relaxation. 

2  Governing  equations. 

2.1  Constitutive  relationship. 

For  every  time  t  £  J  :=  [0,</]  we  consider  the  deformation  of  an  isotropic  viscoelastic  body 
Q,  characteristically  solid,  the  interior  of  which  occupies  the  region  f2  C  TZn ,  where  n  =  2,3, 
with  convex  polygonal  or  polyhedral  boundary  dQ,  and  resisting  the  action  of  a  body  force  /  := 
(/,)"_!  and  a  surface  traction  g  :=  (<Ji)"=i-  The  resulting  displacement  at  a  point  x  :=  (zi)?=i  E 
ft  :=  the  reference  configuration,  is  denoted  by  u  :=  («i)"=1,  and  the  components  of  the 

symmetric  stress  and  strain  tensors  are  denoted,  respectively,  by  cr^  and  £,y ,  where  1  <  i,j  <  n. 
For  the  purposes  of  calculation  these  components  will  be  ordered  in  ^n(n  +  l)-tuplets,  such  that, 
for  n  =  3:- 


a  :=(<rn  C22  033  <’’12  ^2z)T , 

£  :=(£u  £22  £33  £12  £13  £23)T- 


(2.1.1) 

(2.1.2) 
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with  the  reduction  to  lower  space  dimensions  being  obvious.  It  is  to  be  understood  that  (u/t)”=1 
symbolises  an  ordered  n-tuplet,  (uq, . . . ,  w„)T  6  2>",  where  V  is  some  appropriate  set  containing 
the  Wi.  The  strain  tensor  is  obtained  from  the  displacements  via 


_  /  dui  duj  \ 

0  2  +  faTi)  ' 


(2.1.3) 


In  isotropic  linear  elasticity  there  is  a  constitutive  relationship  between  stress  and  strain  given  by 
Hooke’s  law.  In  component  form  this  relation  can  be  expressed  as 

aij(x)  =  A(x)Vt i(x)6ij  +  n(x)eij(u(x)),  (2.1.4)! 

or  equivalently  in  vector-matrix  notation  the  law  can  be  expressed  as 

<r(x)  =  £>(x)e(x),  (2.1.4)2 


where  A  and  fj.  are  the  Lame  coefficients,  <$,j  is  the  Kronecker  delta,  V-  is  the  divergence  operator 
and  D  is  the  constitutive  matrix.  From  (2.1.4)i  and  (2.1.4)2  if  follows  that  the  elements  of  D  are 
given  in  terms  of  A  and  fi.  Thus  for  a  general  compressible  linear  elastic  material,  the  constitutive 
relation  involves  two  independent  scalar  functions  of  position  which,  as  is  given  above,  can  be 
taken  as  A  and  ft.  These  may  be  the  best  pair  of  functions  to  take  for  the  mathematical  neatness  of 
equations  (2.1 .4) i  but  other  functions  are  often  more  useful  in  describing  the  physical  features  of 
a  deformation.  Specifically  we  have  the  bulk  modulus  K,  Young’s  modulus  E,  the  shear  modulus 
G  and  Poisson’s  ratio  v  which  are  given  respectively  by 

k  =  a  +  5'‘'  E  =  q^¥1’  G  =  "/2  “d  "=2(1 Tg}  <215> 

and  for  most  elastic  materials  we  would  expect 

K  >  0,  E  >  0  and  /r  >  0  (2.1.6) 


which  implies  that  —  1  <  i/  <  0.5  but  in  general  we  would  expect  that  0  <  v  <  0.5  which  corresponds 
to  A  >  0. 

In  order  to  arrive  at  a  constitutive  relationship  linking  stress  and  strain  for  a  linear  viscoelastic 
material  at  a  particular  time  t  we  use  hypothesis  1.3.1  in  conjunction  with  a  Hooke’s  law  type 
relation.  We  do  this  as  follows.  We  let  x  6  be  a  fixed  point  and  we  suppose  that  the  strain 
history  at  x  is  of  the  form  shown  in  Fig.  2.1.1.  That  is  we  suppose  that  each  component  of  strain 
is  a  step  function  with  respect  to  the  times  t0  <  t\  <  ...  <  tjv-i-  Specifically  we  have 


f  0,  t<t0, 

\e*(z),  ti<t<ti+i,  i  =  0, 1,  1 


(2.1.7) 


and  thus  the  corresponding  strain  increments  A e{x,ti)  are  given  by 


e°(x), 

e'(x)-£'  :(x), 


i  =  0 

1  <  *  <  N  —  1  . 


(2.1.8) 


The  first  strain  increment  Ae(x,<o)  =  £°(z)  gives  rise  to  a  stress  <7°(x,<)  which  in  the  case  of  a 
viscoelastic  material  decays  the  further  the  time  t  is  from  the  time  t0  when  the  increment  was 
applied.  That  is  mathematically  we  have 


cr°(x,t)  -  D(x,t  -  t0)e(x,t0)}  t>t0, 


(2.1.9) 


where  the  matrix  function  D  in  Hooke’s  law,  which  only  depends  on  position,  is  now  replaced  by 
the  matrix  function  D  which  depends  both  on  position  and  the  elapsed  time  with  the  components 
of  D  decreasing  monotonically  with  the  elapsed  time.  In  (2.1.9),  cr°(x,  t)  is  of  course  also  the  stress 
increment  A<t0(x,  t).  Similarly  the  stress  increment  corresponding  to  the  strain  increment  Ae(x,  t*) 
is  given 


A<t'(x,  t)  —  D(x,  t  -  ti)Ae(x,  <,),  t>ti  . 


(2.1.10) 
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Hypothesis  1.3.1  is  concerned  with  how  these  increments  can  be  combined  to  give  the  actual  stress 
a(x,t)  at  position  x  and  time  t  and  it  states  that  these  increments  can  be  added  in  a  linear  way 
to  give 

AT— 1  N- 1 

<7(a;.  o  =  XI  Aa'(x’ o  =  D(x' 1  ~  <«)Ae(2:’ l') 

«=0  i  =  0 

N-2 

=  Y2  ( D(x,t-ti )  -  D(x,t  —  *i+i))e(;Mi)  +  D(x,t-  tN_ i)e(x,tN_i)  . 
1=0 

(2.1.11) 


Figure  2.1.1  The  step  function  for  the  kl  component  of  strain 


If  D  is  continuously  differentiable  with  respect  to  its  time  argument  then  extending  this  idea 
to  the  limit  as  N  — *  oo  with  — ►  t  and  max,(f,-  -  <i_i)  — ►  0  in  the  usual  way  gives,  using  the  last 
part  of  (2.1.11), 

cr(x,t)  =  D(x,0)e(x,t)  -  f  - }£(x,t)  dr  .  (2.1.12) 

Jto  dT 

In  the  case  of  a  smooth  deformation,  specifically  for  which  de/dr  exists  except  possibly  for  an 
initial  jump  discontinuity  of  e  at  <o,  the  limiting  argument  just  described  gives,  using  the  previous 
part  of  (2.1.11), 

ft 

&{x,t)  —  D(x,t  -  t0)e(x,t0)  +  /  D(x,  t  —  t)  —  (x,  t)  dr  .  (2.1.13) 

Jt0  or 

Note  that  (2.1.13)  could  have  been  formally  derived  from  (2.1.12)  by  integration  by  parts  (and 
vice-versa)  and  also  both  forms  are  implicitly  applicable  only  to  non-ageing  materials  since  the 
stress  relaxation  matrix  D  depends  only  upon  the  relative  or  elapsed  time. 

Next  we  generalise  to  the  case  of  an  ageing  material  by  generalising  the  arguments  of  D. 
Specifically  we  take  D  to  have  the  3  arguments  x,  t  and  r  where  x  is  position,  t  is  the  current  time 
and  t  <  t  is  a  previous  time  and  we  define  the  constitutive  relations 


<r(x,  t)  =  D(x,  t,  0)e(x,  0)  4-  f  D(x,  t,  r)e'(x,  r)  dr 

Jo 

-  D(x,t,t)e{x,t)~  f  D'(x,t,T)e(x,T)  dr, 
Jo 


(2.1.14) 


where  we  have  assumed  for  convenience  that  t0  =  0  and  '  denotes  differentiation  with  respect  to 
t.  In  component  form  the  relation  can  be  written  as 


3-11-1992  16:27  P*.ge  5 


<Tij(x,t)  =  A(i,M)V-u(j,()%+/i(i,(,()£l7(u(i,()) 

A '(x,t,T)V-u(x,T)6ij  +  n'(x,t,T)£ij(u(x,r))  dr,  (2.1.15) 

or 

t)  =  A(x,  t,  0)V'ti(i,  0 )Sij  +  fi(x,  t ,  0)e,-j  (u(x,  0)) 

^(x,t,T)V  u'(x,T)6ij  n(x,t,  r)eij  (u'(x,  r))  dr,  (2.1.16) 

Unless  specifically  stated  otherwise  the  Einstein  summation  convention,  as  used  here,  will  be  used 
in  all  that  follows. 

The  functions  A(x,t,r)  and  fi(x,t,r)  are  the  viscoelastic  analogues  of  the  Lame  coefficients 
that  appear  in  linear  elasticity  (2.1.4).  The  time  dependent  parts  of  these  functions  arise  from  the 
stress  relaxation  functions  characterizing  the  material.  We  must  assume  that  these  functions  exist 
and  that  they,  in  general,  will  vary  for  different  modes  of  deformation,  for  example  bulk  or  shear. 
However,  for  a  non-ageing  material  any  pair  of  these  is,  in  principle,  sufficient  to  determine  the 
time  dependent  parts  of  A  and  (i  via  the  correspondence  principle  [3,4].  The  spatial,  x,  dependence 
of  A  and  p  will  usually  be  of  little  interest  but  is  retained  here  for  generality. 

To  facilitate  the  forthcoming  error  analysis  we  make  the  following,  physically  reasonable, 
assumptions  on  A  and  fi: 

Assumptions  2.1.1. 

i)  For  every  fixed  t  E  Z  and  r  such  that  t  —  t  €  1 

A(x,  t ,  r),  p(x,  t,  r)  £  7im  ^{t  -  rg  Z};  £2(fi)  f)  £TO(fi))  f)  W1'00  |{l-r€  Z};  /^(fi))  , 

where  we  assume  m  may  be  taken  sufficiently  large  for  the  following  analysis  to  remain  valid. 
So  on  a  physical  basis  we  expect  A  and  fi  to  be  well  behaved  in  time,  but  must  allow  for 
material  property  variation  in  space.  This  would  be  the  case  if,  for  example,  two  different 
materials  were  mechanically  joined  in  some  structure. 
it)  For  every  t  E  Z:  0  <  A(x,  t,  r),  /r(x,  t,  r)  Vx  E  fi  and  Vt  —  r  E  Z. 

Hi)  For  every  t  E  Z:  0  <  A'(x,  t,  r),  fi\x,t,r)  Vx  E  fi  and  V<  —  r  E  Z. 

iv)  Causality:  V(  £  I  A(x,  t ,  r),  p(x,  t,  r)  =  0  if  r  >  <.  This  simply  means  that  we  do  not  allow 

future  events  in  the  deformation  history  to  affect  present  behaviour.  ■ 

2.2  Two  weak  formulations:  Pi  and  P2. 

Based  on  the  law  of  conservation  of  momentum  we  consider  the  quasistatic  equations  of  equilibrium 

n  5(7' 

=  /iOM),  *  =  1,-  •  -,ri,  x  E  fi,  tel,  (2.2.1) 

j= l  ax> 

with  the  boundary  conditions 

u(x,  t)  =  0,  x  E  ^  0,  (  £  Z, 

n 

'Y^<rijhj(x,t)  =  flf,(x,<),  x  E  rN,  tel,  i  =  1, . .  .,n 

7  =  1 

Where  TD  [J  VN  =  <9fi  and  n  :=  is  the  unit  outward  normal  to  dCl.  A  weak  formulation  of 

(2.2.1, . . . ,  3)  is  formed  by  taking  the  scalar  product  of  (2.2.1)  with  a  test  function  v  :=  (wj)"=1  £  V, 
where  V<  E  Z 

V  :=  {u  E  (W1^))"  :  t»(x)  =  0  Vx  E  rD|.  (2.2.4) 

Integrating  by  parts  [5]  we  produce  the  weak  formulation  in  which  we  seek  u  E  V  at  each  t  £  Z 
such  that 

/  er(u;  x,  t)  ■  e(u(x))  dx  —  I  fvdx+l  g  vds,V t>  E  V, 

J  n  J  a  Jan 


(2.2.2) 

(2.2.3) 


(2.2.5) 
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where  •  indicates  the  Euclidean  inner  product  and  the  “test  strain” is  given  by 

“>w:=K^+fe)wev'  0£U£n  (2'26) 

Use  of  the  constitutive  relationships  (2.1.15)  and  (2.1.16)  enables  two  alternative  weak  formula¬ 
tions,  Pi  and  P2,  of  the  quasistatic  linear  viscoelastic  stress  analysis  problem  to  be  defined.  For 
notational  convenience,  we  let 

a(A (x,t,r),p(x,t,T),u(x,t,r))  =  I  \(x,t,T)Vu(x,t)V-v{x)  + n(x,t,T)eij{u(x,T))eij(v(x))  dx 

Jn 

(2.2.7) 

which  is  essentially  just  the  usual  bilinear  form  used  to  describe  linear  elasticity.  The  two  problems 
are  then  as  follows. 

Pi:  Find  u(x,<)  £  £2(I;  V)  such  that 

a((A(x,  t,  t),  n(x,  t,  t),  u(x,  t)),  v(x)) 

-  /  a((A/(x,<,r),/j/(x,<,r),u(x,r)),u(x))  dr  =  l(v\t)  Vv  £  V  and  Vt  £  2.  (2.2.8) 

Jo 


P2:  Find  u(x,  t)  £  V)  such  that 

a((A(x,  t,  0),  p(x,  t,  0),  u(x,  0),  v(x)) 

+  f  a((A(x,  t,  r),  n(x,  t,  r),  u'(x,  r)),  v(x))  dr  =  l(v;  t)  Vv  £  V  and  Vt  £  I,  (2.2.9) 
Jo 

where  /(v;  t)  :=  (/(x,  t),  v)n  +  ( g{x ,  <),  v)gn,  and  (•,  )b  is  the  inner  product  on  £2(£).  In  each  case 
we  have  assumed  that  it  is  permissable  to  interchange  the  order  of  spatial  and  temporal  integration. 
In  much  of  the  following  the  x  dependence  of  A  and  p  will  not  be  indicated  explicitly. 

In  order  to  facilitate  the  analysis  we  first  state  some  assumptions,  give  some  notation  and 
state  three  inequalities  which  will  be  required. 

Assumptions  2.2.1. 

i).  For  Pi  and  P2:  f(x,t)  £  £2  (Z;  (£2(fi))n) . 

it).  For  P!  and  P2:  g(x,  t)  £  £2(Z;  {H *  (<9f2))n_1) .  ■ 

Definition  2.2.1.  For  u,  £  7ip(Cl)  we  have 


Definition  2.2.2.  For  v  £  ('Hp(Cl))m  we  have 


Definition  2.2.3.  For  v  €  £2(^;'WP(^))  we  have 


IMI  C,(l,nr) 


a 


\  1/2 

Mlp,n(r)  dr) 
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Lemma  2.2.1.  Korn’s  inequality  [8],  for  v  £  (W1(f2))  ,  and  £(j  defined  by  (2.2.6)  we  have 


£ij(v)tij(v)  dx  >  C||w||i  n, 


for  some  constant  C  >  0  .  ■ 

Lemma  2.2.2.  Continuous  Gronwall  inequality  [9,12],  Let  u,  v,  w  be  piecewise  continuous  non¬ 
negative  functions  defined  on  the  interval  t  £  [0 ,  a] ,  t;  being  non-decreeising.  If,  for  each  t  £  [0,a], 
3 C  >  0,  independent  of  t,  such  that 


then 


u(t)  +  w(t)  <  v 


u(s)  ds, 


u(t)  +  w(t)  <  v(t)  exp(Ct)  . 


m 


Lemma  2.2.3.  Discrete  Gronwall  inequality  [9],  Let  u,  v,  w  be  non-negative  functions  defined  on 

3k  :=  {t  £  H  :  t  =  jk\  j  =  0, 1, . . .,  /;  T  =  Jk} 

with  v(t)  being  non-decreasing  and  k  >  0  being  constant.  If,  for  each  t  6  3k >  3C  >  0,  independent 
of  t,  such  that 


then 


t-k 

u(t )  +  w(t )  <  v(t)  +  Ck  ^  «(t), 

r=0 


u(t)  +  w(t)  <  v(t)exp(Ct)  . 


Remark  2.2.1.  Under  very  mild  restrictions  on  A  and  p  which  are  covered  by  assumption  2.1.1, 
it  can  be  shown  that  a((X,p,  •  ),  •  )  is  a  symmetric  continuous  V-elliptic  bilinear  form  for  each 
(,(-r£l  see,  for  example  [6].  The  V-ellipticity  may  be  shown  by  first  observing  that 

3 C(t)  €  [0, 3]  such  that  \\V-u\\o,n  <  Vu  €  V, 

and  then  using  the  Korn  inequality  (lemma  2.2.1)  to  write: 

I  AV  uV-u  dx  +  /  p£ij(u)eij(u)  dx  >  (X(t)  +£(<)) ||u||i 
J  n  J  n 

where  p (t)  >  0  by  assumption  2.1.1.  ■ 

Assumption  2.2.2.  We  shall  assume  throughout  that  there  exists  a  unique  u(x,t)  €  Ci(l\V) 
that  solves  Pj.  This  in  fact  may  be  shown  using  the  Picard  iteration,  detailed  by  Linz  [13],  and 
exploiting  the  V-ellipticity  of  a((A,  p,  ■  ),  •  ).  ■ 

3  Semi-discrete  formulation. 

3.1  Preliminary  notation. 

In  this  section  we  define  semi-discrete  problems  approximating  Pi  and  P2  by  forming  their  cor¬ 
responding  finite  element  approximations  in  space.  For  this  the  region  fi  is  partitioned  into  Me 
disjoint  finite  elements  such  that  V<  G  I 


Me 

dDh  =dfl  and  tth  :=  (J  Qj1, 

t=l 


(3.1.1) 
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i.e.  the  partition  does  not  change  with  time.  Also  we  define  over  this  partition  the  finite  dimensional 
space 

Vh  :=  peV:  t/V  e  {vp(x))n,  1  <  t  <  Me },  (3.1.2) 

where,  ( [Vp(x ))  denotes  the  set  of  pth  order  complete  polynomials  in  the  variables  (xj )J=1 . 

We  will  approximate  the  weak  solution  u(x,  t)  to  either  of  (2.2.8)  or  (2.2.9)  by  uh(x,  t)  written 
in  terms  of  the  basis  functions  {4>j(x)}jLi  defined  over  .  Specifically  we  have 

Mh 

(uh(x,t)).  =  J2<f>j(x)Lji(uh(x,t)), 

i= l 

where  Mh  >  0  is  an  integer  depending  upon  the  partition  relating  to  the  number  of  nodal  pa¬ 
rameters  (usually  this  will  simply  be  just  the  number  of  nodes.)  and  LJt-  is  the  operator  which 
extracts  the  ith  displacement  component  associated  with  the  jth  basis  function  for  the  finite  element 
approximation.  In  this  way,  for  a  given  i  £l  we  write 

uh(x,t)  =  N(x)U(t),  (3.1.4) 

where  for  any  x  £  Q.  and  any  t  £  1 


( mf  =  Ln(uh(x,t))...LM,n(uh(x,t)))  ennM 


N(x)=[^(x)  ...  *„»(*)]  e  (7f1(u,fi,-))n,nMk. 

<*,(*)  =/*,(*) 

where  I  is  the  identity  on  Hn.  Using  (2.1.3)  we  define  the  approximate  strain  as 

eh(uk(x,t))  =  B(x)U(t) 


(3.1.6) 

(3.1.7) 


(3.1.8) 


where  B(x)  £  (C2(£l))n('n  +  l^2’nM  is  a  matrix  whose  entries  consist  of  derivatives  of  the  basis 
functions. 

3.2  Semi-discretization  of  Pl5  and  error  analysis. 

3.2.1  Semi-discretization. 

The  approximate  stress  is  defined  via  (2.1.15)  and  (3.1.8)  yielding 

crh(uh-,x,t)  =  D(x,t,t)B(x)U(t)—  f  D'(x,t,T)B(x)U(r)  dr.  (3.2.1) 

Jo 

We  then  have  the  semi-discrete  approximation  to  Pi: 

Pi :  Find  uh(x ,  t)  £  C2(X\  Vh)  such  that 

a((A(x,  t,  t),  n(x,  t,  t),  uh(x,  t)),  vh(x )) 

~L  a((X'(x’t'T)’fi'(x’t'T)’uh(x’T^'vh(x^  dr  =  l(vh’t)  £  Vh ,  (3.2.2) 


which  may  be  written  as 


A(t,  t)U(t)  -  [  A\t,  t)U(t)  dr=F(t), 
Jo 


(3.2.3) 


where  A  is  the  positive  definite  (due  to  (2.2.2)  and  assumption  2.1.1)  finite  element  stiffness  matrix 
defined  by 


A(t,s):=  f  Bt (x)D(x,t,  s)B(x)  dx. 
Jn k 


(3.2.4) 
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Remark  3.2.1.  Note  that  (3.2.3)  may  be  written  as 


where 


U(t)  =  G(t)  + 


dr, 


G(t)  :=  (A(t,t))  ' F(t)  £  , 

K(t,r)  =  {A(t,t))~lA'(t,T)  £  (. C2(l))NlN 


for  an  integer  N  depending  upon  the  partition  of  ft  and  the  value  of  n.  The  assumptions  2.1.1  and 
2.2.1  ensure  that  there  exists  U(t)  £  £2(T)  that  solves  (3.2.3),  see  for  example  Tricomi  [10]  where 
the  analysis  for  a  scalar  equation  is  given.  B 


3.2.2  Error  analysis. 

In  this  section  we  derive  a  7f1-error  estimate  for  the  semi-discrete  approximation  uh(x,t),  given 
by  (3.2.2),  to  the  solution  u(x,t)  of  (2.2.8),  for  this  we  employ  the  technique  of  elliptic  projection, 
see  [7].  For  any  t  £  I  the  elliptic  projection  uh(x,t)  £  £i{l>  Vh)  of  u(x,t),  the  solution  to  (2.2.8), 
onto  the  test  space  Vh  is  defined  by 

a((A(x,t,t),fi(x,t,t),u(x,t)-  uh(x,t)),vh)  =  0  Vi>A  £  Vh .  (3.2.5) 

If  we  define 

r)  :=  u  —  uh,  (  :=  uh  —  uh, 

set  v  =  vh  in  (2.2.8),  subtract  (3.2.2),  use  the  projection  (3.2.5)  and  omit  the  implicit  dependence 
on  x  for  clarity  we  have:  ^ 

a((A(M),A‘(M),C(0)>t,h)  =  /  a((A'(^r)»/(<.r).'?(r))>^)  dr 

Jo 

+  /  a((x'{t,T),fi'(t,T)X(r)),vh)  dr.  (3.2.6) 

Jo 

Since  we  seek  to  estimate  the  error  ||u  —  using  the  above  definitions  we  have  immediately 

that,  for  t  £  1, 

||u(:M)-uft0M)l|i,o  <  lb(M)lkn  +  ||<OM)lkn,  (3.2.7) 

and  the  task  is  now  that  of  bounding  the  terms  r?  and  (  in  the  norm  ||  •  ||i  n- 

Now,  for  any  t  £  J,  we  choose  vh  =  C(t)  in  (3.2.6)  and  observe  that  remark  2.2.1  gives  a 
C\(t)  >  0,  thus  allowing  (3.2.6)  to  be  written  as 

C'i||C(0lli,n  <  f  “((^d)./<'(',f),lW),((())  dr 
Jo 

+  /  a((A'(£r)>*1'(£r)>C(7')),<(<))  dr-  (3.2.8) 

Jo 

We  proceed  by  noting  the  following 
Lemma  3.2.1.  Define  the  symbol  6(-  ,  •  ,  •  )  by 

K4>,  V-T7,  V-0  ■=  f  f  4>(x,  <>  «)V  '7(z,  s)V-C(x,  t)  dxds 
Jo  Jn 


where 


<j>(x,t,s)  £  £2((0,<);£oo)  f°r  ea°h  t  e  ^ 

r?(x,s)££2((0,0;(7i1(fi))") 

C(x,<)  £  (Tf^ft))"  for  each  t  £  X. 
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Then 


b(4>,  V  rj,  V-C)  <  ||^(s)||Q,(o,t)|K(Olli,n||f?(s)||£2((o,i);(w-(«))n) 


Vt  £1, 


where 


<^(s)  =  ess  sup  if>(x,  s ). 

rgn 


Proof.  Interchanging  the  order  of  integration  and  liberally  applying  the  Cauchy-Schwartz  inequal¬ 
ity  gives 


b(<t> 


,  v-77,  V-C)  =  J  (J  <Hs)V-n{s)  ds)v<(0  dx 


2  -I  1/2 

d>(s)V-r/(s)  ds  )  dx 


T  1/2 


<IIC||i,n[|nll^,(0,()  l|V  ^||2  (0  t)  dx 

<  Iloilo, (0,t)IKI|i>a||»?||i:j((0,i);(w>(n))“)-  ■ 

Noting  the  definition  of  a((A,  //,-),-)  and  applying  lemma  3.2.1,  coupled  with  an  almost  identical 
argument  for  the  e  terms  (i.e.  consider  b(<j>,  e{r?),  £(())  in  the  lemma),  we  write  (3.2.8)  as 


||C(0llx,n  <  +  ||C(i‘)lk2((o,t);(w>(«))»))  (3.2.9) 


with 


C2  —  Cxl  ^||A'||0](o,t)  +  9||Allo,(o,t)) 


Cx  :=  inf  Ci. 

t-re[0,t] 


Squaring  both  sides  of  (3.2.9)  and  using  the  inequality 

2afe  <  a2  +  62  Va,  b  £  7£, 
we  have  by  Gronwall’s  inequality  (lemma  2.2.2) 


(3.2.10) 

(3.2.11) 


(3.2.12) 


HC(0lli,n  <  2C2||r7(r)|||2((01).(Ki(n))„)  +  2C2||C(i’)||ic2((0,t);(«i(n))«) 

IIC(0lli,n  <  G3||r?(r)||£:2((o,t);cw'(n))")exp(C22<),  where  C3  :=  y/2C2.  (3.2.13) 


Assuming  a  suitable  discretisation  of  the  domain  and  using  assumptions  2.2.1  we  have,  see 
[6],  constants  >  0  and  a  >  0,  depending  upon  the  regularity  of  the  weak  solution,  the  degree  p 
of  the  approximating  polynomials  and  the  mesh  partition  (i.e.  the  space  Vh .),  such  that 


where 


IMs.Olkn  <  C  inf  ||u  -  t/*||i,n 

<  C,fia|u(x,0li+a,n  Vt£l, 


h  max  diam  flj1 . 
l<i<M,  1 


(3.2.14) 


Hence,  we  have  the  following  semi-discrete  error  estimate: 


Theorem  3.2.1.  Assuming  that  u  £  C2(X\  V)  {^\C2(j]  (?fI+“(0))n)  (\Coo(l\  (?f  1+a(f2))")  and 
that  uh  £  Cx(X\Vh)[>\C2{T\Vh),  and  observing  assumptions  2.1.1,  then  for  the  error  associated 
with  P  j ,  the  semi-discrete  aproximation  to  Pi,  Vt  £  I  B  C  >  0  depending  upon  u  but  independent 
of  h  such  that 


||u(x,<)  -  uh(x,  <)||i,n  <  Cha. 


3-11-1992  16:27  Pa-ge  11 


Proof.  Use  (3.2.7)  with  (3.2.13,14)  to  write  at  1  6  I: 

||«(x,<)-«h(*iOlli.n  <Cr,/ia|u(x,<)|i+a,n+C'3C^/i°|w(a:,T)|jC2((0  t);(WH<»(n))->)exp(C2t/2) 

<C,/i"|«(x,  r)|£oo(i;(«i+«(n))»)  +  CsCfl^K*.  r)U2(X;(«i+«(n)).)  exp(Cf</2) 
<Cha 


where 


and 


C-i  :=  sup  Ci. 

t-rel 


C3  :=  sup  C3, 

t-rez 


Tie  existence  of  C?  a nd  C3  are  guaranteed  by  assumption  2.1.1  and  remark  2.2.1. 


3.3  Semi-discretization  of  P2,  and  error  analysis. 

3.3.1  Semi-discretization. 

Again  using  (3.1.8)  but  this  time  with  (2.1.16)  we  define  the  approximate  stress  through 

ah(uh;x,t)  =  D(x,t,0)B(z)U(0)+  I  D(x,t,r)B(x)U\r)  dr.  (3.3.1) 

Jo 

In  which  case  we  have 

Pj:  Find  uh(x,  t )  6  Vh  for  each  t  E  1  such  that 
a((A(x,  t,  0),  p(x,t,  0),  uh(x,  0)),  vh(x)) 

+  f  a((\(x,t,T),p(x,t,r),(v.h(x,r))'),v\x))  dr  =  l(vh-,t)  Vuh  €  (3.3.2) 

Jo 

which  may  be  written  as 

[  A(t,r)U'(r)  dr  =  b(t),  (3.3.3) 

Jo 

where 

b(t):=F(t)-A(t,0)U(0).  (3.3.4) 

The  initial  condition  f/(0)  is  found  via  the  equation  6(0)  =  0,  i.e.  by  solving  a  problem  of  linear 
elasticity  at  t  =  0.  In  this  case  we  have  produced  a  Volterra  integral  equation  of  the  first  kind  for 
U'(t). 

3.3.2  Error  analysis. 

As  a  rule  Volterra  integral  equations  of  the  first  kind  are  rather  more  difficult  to  deal  with  than 
those  of  the  second  kind  (see  for  example  [10,13]),  moreover  once  the  initial  condition  has  been 
determined  (2.2.9)  is  actually  an  integrodifFerential  equation  for  u(x,<)  of  a  non-standard  type. 
We  have  not  attempted  a  direct  semi-discrete  error  analysis  for  P2  here,  but  choose  instead  to  rely 
on  the  observation  that  if  the  numerical  algorithm  promised  by  P2  is  to  be  useful  then  u'(x,t )  is 
required  to  exist,  at  least  in  the  weak  sense,  and  so  we  may  simply  state: 

Theorem  3.3.1.  If  the  conditions  required  for  theorem  3.2.1  are  satisfied  with  the  further  con¬ 
dition  that  u(x,t)  €  W1'1(I;y)nT2(l;V)  Cl^2(l;(n1+a(Q))n)  f]£00(l-,(n1+a(Q))rl)  -a  suffi¬ 
cient  condition  for  this  extra  regularity  in  time  may  be  provided  by  strengthening  assumption  2.2.1 
to  f(x,  t )  G  C1  (J;  (£2(fl))")  a  nd  g(x,  t)  G  Cl(l ;  (H  2(<9ft))n_1),  then  for  the  error  associated  with 
P2,  the  semidiscrete  approximation  to  P2,  3 C  >0  Vi  G  I  such  that 

||«0M)~  uh0M)||i,n  <  Cha, 
and  this  C  is  the  same  as  that  appearing  in  theorem  3.2.1. 

Proof.  Integrate  (2.2.9)  by  parts  to  yield  (2.2.8)  and  (3.3.2)  to  yield  (3.2.2),  because  of  the 
uniqueness,  assumption  2.2.2,  the  result  of  theorem  3.2.1  may  be  applied.  ■ 
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4  Fully  discrete  formulation. 

4.1  Preamble. 

To  provide  fully  discrete  approximations  to  Pj*  and  P!j  we  use  respectively  the  trapezoidal  rule 
for  numerical  integration  and  a  type  of  product  integration  rule  coupled  with  a  finite  difference 
replacement.  In  the  ensuing  analysis  we  assume,  in  both  cases,  a  constant  time  step  k  :=  t,  — 
fj_i,  i  =  1, . . . ,  I-  The  time  domain  1  is  discretized  into 

lk  :=  {0  =  t0,  ■  •  •  ,<i,  •  •  ■  ,</  €  H  :  and  t,  >  h_i  for  1  <i</}  Cl.  (4.1.1) 

4.2  Full  discretization. 

4.2.1  Full  discretization  of  Pj. 

In  order  to  present  a  fully  discrete  scheme  for  Pi  we  need  to  discretize  P*  in  time,  for  this  we 
utilise  the  trapezoidal  rule  for  numerical  integration,  the  generic  form  of  which  is: 


r tj  _i s  fc. 

8(k,tj)=  y(x)  dx-2_J-A(yi-i+yi), 

JO  i=l  1 

(4.2.1) 

where 

y<  ■=  y(xi), 

ki  ■—  Xi  —  Zi_i,  1  <  i  <  j  <  I , 

(4.2.2) 

Then  by  setting 

k  =  max  {kA, 
i<i<j  1  ’ 

(4.2.3) 

and  assuming  appropriate  continuity  of  y,  8  may  be  estimated  by 

\8{k,tj)\  <  Cjk2\y”(Z)\,  for  some  ^  £  (0, <j). 

(4.2.4) 

Thus,  we  replace  the  integral  occuring  in  (3.2.3)  with  the  trapezoidal  approximation,  neglect 
the  error  term  8,  and  arrive  at: 

Pk,k:  Find  ( uh(x ))^.  6  for  each  tj  £  Xk  such  that 


+  =  (4.2.5) 

?= i 

Rearranging  this  equation  to 

j  - 1  . 

,  tq)Uq  +  A'(tj ,  tq~i)Uq-i 

and  taking  kj  small  enough  for  the  left  hand  side  term  to  be  invertible  allows  successive  solution 
for  Uo,U\,  ...  ,[//,  thus  yielding  the  approximation  to  uh(x,t)  by  substituting  Uj  for  U(tj)  in 
(3.1.4)  giving  (uh(x))j .  We  note  here  that  the  trapezoidal  method  does  not  require  starting  values 
and  is  not  restricted  to  constant  k  as  higher  order  methods  would  be.  Also,  since  the  scheme  has 
a  repetition  factor  of  1  it  is  numerically  stable  in  sense  of  Linz,  see  for  example  [13,14]. 

4.2.2  Completion  of  the  error  estimate  for  Pj. 

To  build  upon  the  semi-discrete  error  estimate  of  theorem  3.2.1  we  will  denote  the  fully  discrete 
solution  to  Pj’fc  at  each  <t  £  lk  by  ( uh(x) ).  and  also  define  d(x,  tr)  :=  uh(x,  tr)~  ( uh(x))r  V<r  £  Xk. 
We  then  have 

llU0Mr)  -  («k(*))r||  <  IK»,<r)  -  u'*(z,tr)||  +  P0Mr)||.  (4.2.6) 

As  the  first  term  on  the  right  hand  side  of  (4.2.6)  has  been  bounded  in  theorem  3.2.1  we  now 
seek  an  estimate  for  9  in  the  ||  !]i,u  norm.  Assuming  a  constant  time  step  k,  (4.2.5)  may  be  written 


-  —A'itjJj^Uj  =  F,-  +  y 
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as: 

Find  ( uh(x))r  £  V*  for  each  tr  £  1k  such  that  6  V* 

r 

a((X(tT,tr),fi(tr,tr),(uh)r),vh )  -  k'^wrqa((X,(trJq),n'(tr,tq),(uh)ll),vh)  =  Fr,  (4.2.7) 

1=0 

where  the  wrq  are  the  weights  associated  with  the  trapezoidal  rule  (4.2.1).  Subtraction  of  (4.2.7) 
from  (3.2.2)  evaluated  at  t  =  tr  gives  (using  (4.2.1)) 

r 

a({X(tr,tr),n(tr,tr),0(tr)),vh )  -  6(k,tr)  +  k^2wrq  a((A'(fr,  tq),  /i'(tr,  tq),  9(tq)),  vh)  (4.2.8) 

9=0 

and,  assuming  certain  smoothness  in  time,  specifically  uh(x,)  £  C2(l)( for  this  it  is  sufficient  to 
strengthen  assumption  2.2.1  such  that  l(v;  ■)  £  C2(l),  see  for  example  [13].) 

6(k,tr)  =  Crk2-^a((X'(tr,f),fi'(tr,f),uh(x,T)),vh),  t  £  [0,fr],  (4.2.9) 

We  choose  vh  =  9r  9(tr)  £  Vh  and  again  use  the  V-ellipticity  (remark  2.2.1)  on  the  term  on  the 
left  hand  side  of  (4.2.8)  giving  a  constant  Ci(f)  >  0  such  that 


CiOr)||0r Hi^n  <  &{k,tT)  +  k  'y  ]  wrq  a((X  (tr ,tq),  [i  (tr,tq),  6 q),  0r) .  (4.2.10) 

i= o 

To  bound  the  right  hand  side  of  (4.2.10)  we  observe  that 

[  X \x,tr,tq)V-9(x,tq)V-6{x,tr)  dx  <  f  |A'(x, <r, f,)|  |V-0(x,t,)|  |V-0(x,tr)|  dx 
J  n  J  n 


<  Cqr( X)  110,11!, n  ||0r||i,n, 


(4.2.11) 


and 


where 


[  /A x,tr,tq)eij(9q)eij{9r )  dx  <  f  \^'(x,tr,tq)\  |e,;(0,)|  |e,;(0r)|  dx 
J  n  Jn 

<  9 C,r(/i)  ||0,||i,n  ||<M|i,n-  (4.2.12) 


0  <  C,r(A)  ||A'(x,  tr,  i,)||o,oo,n, 
0  <  C,r(/i)  :=  \\n'(x,  tr,  <,)||o,oo,n- 


The  strict  inequality  is  a  consequence  of  assumption  2.1.1.  Then  combining  (4.2.10-12)  gives 

r 

Ci(*r)||0r||?,n  <  l*(Mr)|  +  k\\9r || , ,n  £  wrq (Cqr (. X)  +  9 C,r(/i))  ||fl, 111,. n 

5  =  0 


with 


<  |6(fc,fr)|  +  C2(fr)^||0r||l,n  ^2 

?  =  0 


C2(tr)  ■=  ess  sup 


{hr? |  |C,r(A)  +  9C,r(//)|  j. 


(4.2.13) 

(4.2.14) 


Turning  to  the  quadrature  error,  6,  we  see  from  (4.2.9)  that 

m,tr)\  <  |Cr|*2(|  JjX'(tr,T)Vuh(T))"V-9r  dx\+\jjii'(triT)elJ(uh(f)))"e,1(0r)  dx 


<  \Cr\k2  ^||(A'(<r,  f)V-uh(r))  ||0,n  +  \\(n'(tr,f)eij  (u,*(f))),,||oin  )  ||*r||i,n 

<  C3(tr,  uh)fc2||0r||i,n. 


(4.2.15) 
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Using  this  in  (4.2.13)  we  have 


Ci (tr)||0P||i.n  <  C3(tr,uh)k2  +  C2(tr)k  IKIli.n 

?=o 

or,  for  k  small  enough,  by  rearrangement  this  is 


r  — 1 


(Ci(<r)  -  *C2(<r))||*r||1>n  <  c3(tr,  uh)k2  +  Ca(tr)t  ll*,lli.n, 

9=0 


therefore 


||0r||i,n<  C5(«fc)Jb2  4- C4Jb5ZH^||1,n  Vfr  e  Ik, 

q=0 


(4.2.16) 


where 


ess  sup 
tei 


C3(t,uk) 

Ci(t)-kca(ty 


We  may  now  bring  these  results  together  and  state  the  following  theorem  for  the  discrete  problem 

PM; 


Theorem  4.2.1.  If  the  conditions  of  theorem  3.2.1  are  satisfied  together  with  the  additional  re¬ 
quirements:  i)  k  <  Ci(t)/C2(t)  Vi  £  X,  ii)  assumption  2.2.1  is  strengthened  so  that  g(x,  •),  f(x,  ■)  £ 
C2(  1).  Then  for  the  error  associated  with  Pj'k,  the  fully  discrete  approximation  to  Vtr  £ 
Xk  3C  >  0  depending  upon  u(x,t)  and  uh(x,t)  but  independent  of  h  and  k  such  that 

IK*.  <r)  -  (uh(2:))r||i,n  <  C(ha  +  k2)  Vfr  £  Jk. 

Proof.  Apply  Gronwall’s  inequality  (lemma  2.2.3)  to  (4.2.16)  to  give 

ll^r||i,n  <  C3k2  exp(Cqti)  Vtr£lk, 


where  the '  denotes  suprema  over  1  of  the  appropriate  quantity,  and  then  combine  this  result  with 
theorem  3.2.1  using  (4.2.6)  to  obtain  the  desired  result  m 


3  Full  discretization  of  P2 


'  4.3.1  JEuULdiscretization. 

r'(  Witl\(4.2^2,3))still  in  force  we  obtain  a  fully  discrete  analogue  of  (3.3.3)  by  employing  the  following 

•  finite  difference  replacement  [11]: 


where 

A Ui  Ui  -  Ui-\,  1  <  t  <  I,  with  Ui(x)  ss  u(x,  ti). 
Applying  this  to  (3.3.3)  we  obtain: 


(4.3.1) 

(4.3.2) 


P^’*:  Find  (uh(x)).  £  Vh  for  each  tj  £  lk  such  that 

J  Af/  ftq 

Yl-iT1  A(tj,r)  dr  =bj,  1  <  j  <  /,  (4.3.3) 

q= l  JU-i 

with 

Uo  =  (A(0,0))-1Fo,  (4.3.4) 

where  (4.3.4)  follows  from  (3.3.4).  Rearranging  (4.3.3)  U,  is  given  explicitly  in  terms  of  known 
quantities.  Note  that  this  method  of  discretization  corresponds  closely  to  the  midpoint  rule  for 
Volterra  first  kind  integral  equations  which  is  known  to  have  desirable  convergence  and  stabilty 
properties  [13].  The  numerical  results  presented  in  §5  bear  this  observation  out.  Again  we  note 
that  the  method  is  not  restricted  to  constant  k  as  higher  order  methods,  obtained  by  higher  order 
finite  difference  approximations,  would  be. 
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4.3.2  Completion  of  the  error  estimate  for  P2. 

We  again  use  9  as  defined  in  §4.2.2  and  employ  a  constant  time  step  k  to  provide  a  bound  for  6r 
in  the  norm  ||  •  ||i,n-  We  commence  by  evaluating  (3.3.2)  at  t  =  tr  £  lk  and  subtract  (4.3.3)  from 
this,  giving 


/  a((Htr,r),n(tr,T),xq),vh)  dr 
<J=1  Jtq-l 


=  0  Vi/1  e  vh 


where 


Xg  ■=  (uh(X,T))'  -  T  €  (tq-\,tq) 

A Uq  :=  Uq  —  Uq—i . 


(4.3.5) 


(4.3.6) 

(4.3.7) 


Now,  assuming  sufficient  smoothness  with  respect  to  time,  we  use  Taylor’s  series  about  the  mid¬ 
point  of  each  time  interval  (<?_i,<?).  Specifically  we  assume  ti(x,  •)  6  C3(l)  (i.e.  we  strengthen 
assumption  2.2.1  to:  f(x,  •),  g(x,  •)  €  C3(I).)  and  define  for  1<9<  r  and  r  6  (tq-\ , tq) 


(u  (tq)  u  (tq- 1)) 


%,  r)  =  (uh(r))'  -  — JlA-Hi-il 

and  note  that  Taylor’s  series  with  the  Lagrange  form  of  the  remainder  gives 


(4.3.8) 


6{q,r) 


al.»(T) 


V"  ,  ( T~Cq )2 


=  (r  -  Cq)(uh(Cq))"  +  ^-^(«A(ai ,,(r)))'"  -  (4.3.9) 

=  (tq-  1  +  tq)/ 2 

=  cq  -f-  w(r  —  cq)  for  some  u>  £  (0, 1) 


012, q  €  (fj-1,  <?). 

Hence,  since  uh(x,tq )  =  t/?(x)  +  6q  the  function  is  given  by 


Xq  == 


_ Qq  ~  Qq-  1 


+  5(q,r)  for  t  £  (tq_l}tq),  1  <  q  < 


(4.3.10) 


Now  we  choose  vh  —  (9(x))r  in  (4.3.5)  giving 

r  rtq 

YZ  /  a((Ktr,r),n(tr,T),(9q  -  9q-i)/k  +  6(q,T)),er)  dr  =  0, 

5=1 

expanding  the  definition  of  a(( A,/z,  •  ),  •  )  given  by  (2.2.7)  we  write  this  as 

tV  /  V-(5,  -  0?_i)V-^r  /  A  (tr,r)  dr  +  £ij(9q  -  0q-i)eij(9r)  f  n(tr,r)  dr  dx 

q  =  1  /‘l-I 

+  YZ  f  a((Ktr,T),^(tr,T),6(q,r)),er)  dr  =  0.  (4.3.11) 

<?=i 

Rearranging  the  sum  gives,  equivalently 

j;  /  V-9rV-9r  f  X(tr,T)  dr  +  £ij(9r)£ij(9r)  [  pi{tr,r)  dr  dx 
*■  J  n  d  t  r 1  /fr-l 

+  -  V]  f  V-(?rV-0,  I  A (<r,r)dr-  f  A(<r,r)  dr  dx 

k  q=lJn  I  •/«,  J 

+  tYZ  f  £ij(er)£ij(0q)  f  fi(tr,T)dr-  [  n(tr,r)dr  dx 

q—  1  Jtq-1  dtq 

r 

+  'll  a((Ktr,T),v(tr,T),6(q,T)),6r)  dr  =  0.  (4.3.12) 

1  •'t,-! 
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Due  to  the  fact  that  A  and  p  may  be  assumed  continuous  in  time  (assumptions  2.1.1)  this  may  be 
rewritten,  by  liberally  applying  the  mean  value  theorem,  as 

j  A(tr,Tir)  V-#rV-dr  +  j  T2r)£tj(@r)£ij(&r)  dx 

Jn 

r~ 1  t 

-  CkkY]  /  A,(tr,rlr,)V-0rV-0?+/(<r,f2r?)£.;(^r)^(0?)  dx 

q=lJn 

+  ^2  j  a(mtr,T),Ktr,r),6(q,T)),er)  dr  =  0,  (4.3.13) 

9  =  1  dtq- l 

where  0  <  C*  <  2.  Now,  A(tr,fjr),  fi(tr, f2r)  >  0  by  assumption  2.1.1,  so  observing  remark  2.2.1 
3Ci(tr)  >  0  such  that 


C,l(tr)||^r||i  n  —  Ckk  'y  /  A  (<r,  fir?)V-0r  V  dg  +  fJ,  (<r,  T2rq)Cij(9T)eij(6ij) 

t=iJn 

f  «((A (tr,T),tx(tr,T),6(q,T)),6r)  dr. 


dx 


/-j  / 

9=i  ■'*»- 


(4.3.14) 


Taking  absolute  values  we  express  (4.3.14)  as 

r— 1 


C1(tr)||*r||?,n  <  E  P*IMN|i,n(|A'(tr,  hr,)\  +  9|/i'(tr,  f2r?)|) 


1- 1 


(4.3.15) 


+  E  /  a((A(<r,r),/i(tr,r),%,r)),er)  dr  . 

7=i 

Turning  now  to  the  term  involving  the  truncation  error  S.  We  note  that  from  (4.3.9), 

[  a((A(<r,r),/%r,r),%,r)),0r)  dr  =  f  a((A(fr,  r), n{tr,  r),  %,  r)),  9r)  dr  (4.3.16) 

Jtq-\  Jtq- 1 


where 


%,  r)  :=  („‘(ftM(r)))'"  -  g  (ti*(a2.,))' 


(4.3.17) 


(4.3.18) 


as  a  consequence  of  the  elementary  result 

f  (r  —  c9)dr  =  0 

Thus  we  have  for  each  component  of  the  sum  in  (4.3.15) 
a((A(<r ,  r),  p(<r,  t),  %,  r)),  0r)  dr 

<  /  |V-0r|  /  | A(ir ,  r)|  |V-6(q,  r)|  drdx  +  f  |e,-y(0r)|  /  !M<r,r)|  |e<y  (%>  r))l  drd;c 

J  CL  Jtq—i  J  Cl  Jtq  —  i 

fc|A(tr,r?)|  f  |V-<9r|  |V-%,f,)|  dx  +  k\fi{tr,  r?)|  /  |eij(^r)|  \etj  (%,  r,))|  dx 
Jn  Jn 


< 

<  *:|a(4,t,)|  ||%|i,n||%%||i,a  + %(*r%)|  ||0r||i,n||%%)||i,n 
From  (4.3.17)  we  have  the  bound 

11%,  *)||i,n,  ||%,  r)||i,n  <  Cjk 2  Vr£l# 


(4.3.19) 


(4.3.20) 
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Now  if  we  take  appropriate  suprema  over  the  interval  [to,tr]  we  may  write  (4.3.15)  as 

r  —  1 

IKIkn<C2^  +  C3^Ell^!k«  VtreX*.  (4.3.21) 

9=1 

Hence  we  have  the  following 

Theorem  4.3.1.  If  the  conditions  required  by  theorem  3.3.1  are  satisfied  and  further  if  assumption 

2.2.1  is  strengthened  to  f(x,  •),  g(x,  •)  €  C3(l)  then  for  the  error  associated  with  P^'*,  the  fully 
discrete  approximation  to  P2,  Vtr  6  lk  3 C  >  0  depending  upon  uh(x,  t)  but  not  on  h  and  k  such 
that 

IK*,  tr)  -  («'>(x))ri|i,n  <  C(h°  +  k2)  Vtr  €  Ik- 
Proof.  Apply  lemma  2.2.3  to  (4.3.21)  to  obtain 

Pr||i, n  <  C2tik 2  exp(C3</)  V<r  £  lk, 

and  use  (4.2.6)  in  the  context  of  Pij’*.  ■ 

5  Numerical  experiments. 

5.1  Preamble. 

In  this  final  section  we  predicate  the  theoretical  convergence  rates  given  in  theorems  4.2.1  and  4.3.1 
by  computing  the  approximate  solution  to  simple  model  problems  for  which  the  analytical  solutions 
are  known.  These  solutions  are  found  via  the  correspondence  principle  [4],  which  enables  the 
viscoelastic  solution  to  be  derived  from  the  corresponding  linear  elastic  solution  by  using  integral 
transform  techniques.  The  correspondence  principle  is  only  applicable  to  non-ageing  materials. 
Accordingly  we  set,  neglecting  the  x  dependence, 

A(t,r)  =  A(<  -  r),  (5.1.1) 

p(t,r)  =  p(t  -  t).  (5.1.2) 

We  shall  assume  that  the  time  dependent  part  of  A  and  p  may  be  modelled  adequately  by  expres¬ 
sions  of  the  form: 

N 

^(<)  =  ^Cie"0,<.  (5.1.3) 

i  =  0 

which  are  suggested  by  the  appearance  of  stress  relaxation  curves  for  real  viscoelastic  materials, 
see  the  data  for  “Maranyl”[2]  in  §5.4  for  example,  or  the  curves  given  in  [1],  This  form  for  A  and  p 
introduces  an  important  simplification  into  the  numerical  algorithms  (4.2.5)  and  (4.3.3).  Usually 
in  solving  these  kinds  of  integral  equation  problems  the  solutions  obtained  at  each  time  step  up 
until  the  current  one  must  be  stored  in  order  to  compute  the  history  integral.  This  can  be  very 
demanding  on  computational  storage.  However,  for  history  kernels  of  the  form  (5.1.3)  and  for  Pj'* 
we  exploit  the  relation: 

f  '  e-a,(t,-T^(r)  dr  =  f  dr  +  [  '  e-a'(-t'~T) g(r)  dr  (5.1.4) 

Jo  Jo  Jt,-i 

for  any  function  g,  i.e.  for  each  term  of  the  sum  in  (5.1.3)  the  history  integral  at  time  tj  may 
be  obtained  simply  by  attenuating  the  result  at  the  previous  time,  tj- 1,  by  an  amount  dependent 
upon  the  current  time  step.  Thus  only  N  +  1  solution  vectors,  per  viscoelastic  function,  need  be 
stored  in  order  to  evaluate  the  history  term  completely,  see  [11],  A  similar  observation  can  be 
made  for 

In  the  final  subsection  we  compare  the  numerically  computed  results  with  the  experimental 
stress  relaxation  and  creep  data  for  a  real  engineering  material  [2],  which  we  assume  to  be  non¬ 
ageing.  Also,  following  the  manufacturer’s  recommendations  we  take  A (<)  oc  p(t),  the  constant  of 
proportionality  being  obtained  from  the  linear  elastic  definitions  of  the  Lame  coefficients,  which 
are  given  usually  in  terms  of  Young’s  modulus  and  the  Poisson  ratio,  the  latter  of  which  will 
be  constant  when  this  proportionality  exists  between  the  viscoelastic  functions  and  the  former  is 
embodied  in  the  relaxation  functions  themselves. 


3-11-1992  16:27  P*g«  18 


5.2  Benchmark  analytical  solutions. 

5.2.1  Linear  elastic  solutions. 

For  the  problem  (2.2.1,. .  .,3)  we  take  n  —  2  and  define  the  domain  Q  to  be  the  rectangle 

:=  {  x,y  £1Z  :  0  <  x  <  /,  0  <  y  <  2k,  l,  k  >  0 } ,  (5.2.1) 

and  present  two  linear  elastic  solutions  to  the  problems  arising  when  Q  has  certain  body  forces 
and  surface  tractions  imposed  over  it  and  on  it  (resp.).  (Although  l  and  k  have  been  used  in 
previous  sections  with  quite  different  meanings,  the  usage  in  (5.2.1)  is  essentially  confined  to  this 
section  (5.2),  definite  values  of  /  and  k,  in  this  sense,  being  assigned  in  all  the  following  work.) 
The  first  of  these  solutions  has  displacements  that  vary  quadratically  in  x  and  y,  so  that  if  the 
finite  element  employed  in  the  computation  is  composed  of  quadratic  basis  functions  the  spatial 
behavior  of  the  displacements  may  be  captured  exactly  by  the  numerical  scheme  thus  isolating 
the  time  discretization  error.  The  second  solution  is  completely  artificial  in  that  the  body  forces 
are  dependent  on  the  Poisson  ratio  of  the  material,  for  our  purposes  this  is  immaterial  since  we 
are  interested  only  in  the  behaviour  of  the  numerical  solution,  the  displacements  arising  in  this 
case  are  sinusoidal.  The  domain  Q  is  shown  schematically  in  figure  5.2.1,  where,  in  the  previous 
notation 


(z,y)  :=  (ii,x2),  (5.2.2) 

(ue.ve)  :=  («i,u2),  (5.2.3) 

(X,Y)  :=  (/i,/2).  (5.2.4) 


The  subscript  E  denotes  linear  elastic  solutions.  When  this  subscript  is  dropped  in  §§5.2. 2, 3  the 
solutions  will  be  for  the  corresponding  linear  viscoelastic  solution. 


Figure  5.2.1 


i)  Quadratic  solution. 

To  obtain  a  quadratic  variation  of  displacement  we  impose  upon  Q  the  following  boundary 
conditions 

ke(0,  k )  =  vjs(0,  k )  =  ve{1,  k)  =  0,  (5.2.5) 

and  assume  that  a  spatially  constant  body  force,  X(t),  acts  whilst  Y(t)  =  0.  Then  we  have  the 
single  non-zero  surface  traction: 


<7i  =  lX(t)  V(x,y)  G  {x  =  0,  0  <  y  <  2k). 


(5.2.6) 
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Under  these  conditions  the  solution  to  (2.2.1)  under  the  assumption  of  plane  stress  is  given  by 


UE(x,y,t)=  2^(2/a:  z2)  v(y  k)2|- 

(5.2.7) 

ve(z,V,  t)  =  g\x  l)(y  k), 

(5.2.8) 

<TU  =  X(t)  (/  -  x), 

(5.2.9) 

<X  1 2  =  ^22  =  0  • 

(5.2.10) 

Rotating  Q  through  90°  clockwise  we  see  that  this  problem  corresponds  to  the  physical  situation 
of  a  plate  suspended  from  the  midpoint  of  its  uppermost  edge  and  deflecting  vertically  under  the 
action  of  its  own  weight.  In  these  expressions  E  is  the  Young’s  modulus  of  the  material  and  t/  <  1/2 
the  Poisson  ratio. 


ii)  Trigonometric  solution. 

As  remarked  above,  this  is  an  artificial  solution  because  the  body  forces  depend  upon  the 


Poisson  ratio  of  the  material.  We  set 

X(t)  =  cos  <j>x  +  £2xL(t)sin£(y  -  1/2),  (5.2.11) 

Y(t)=  cos£(y-  1/2  )  +  <f>2(y-  l/2)N(t)  sin<t>x,  (5.2.12) 

and  employ  the  essential  boundary  condition  (5.2.5).  Under  the  assumption  of  plane  stress  it  can 
be  shown  that  the  solution  of  (2.2.1)  is  now  given  by 

UE(x,y,t)  =  x sin  £(y  -  1/2),  (5.2.13) 

VE{x,y,t)  =  2(1  +^)jV(<)(y  _  1/2) sin <f>x,  (5.2.14) 

2 

ern(x,y,t)  =  _  ^  L(t)sin£(y  -  1/2)  +  vN(t)s'm<j)x  ,  (5.2.15) 

o~22 (x,y,t)  =  N(t)  sin  <f>x  +  vL(t)  sin  £{y  -  1/2)  ,  (5.2.16) 

crn(x,  y,t)  =  £xL(t)  cos^(y  -  1/2)  +  <j>{y  -  1/2)1V(<)  cos <f>x,  (5.2.17) 

when  the  tractions  are  applied  consistent  with  these  expressions  for  <t,j. 

5.2.2  Linear  viscoelastic  solutions  (non-ageing). 


Linear  viscoelastic  solutions  for  non-ageing  materials  may  be  generated  from  particular  types  of 
linear  elastic  solutions  via  the  correspondence  principle  [3,4].  To  give  a  simple  example  of  how  this 
may  be  achieved  we  firstly  assume  that  the  time  dependent  parts  of  A  and  y  are  proportional,  i.e. 

A(<)  oc  y(t),  (5.2.18) 


this  is  tantamount  to  assuming  a  constant  Poisson  ratio.  In  terms  of  the  D  matrix  of  (2.1.12)  or 
(2.1.13)  this  means  that  we  have,  omitting  again  the  dependence  on  x,  the  representation 

D(t)  =  D0<p(t),  (5.2.19) 

where  Do  is  a  constant  matrix.  Thus  we  have  the  relation 


D0d(t  -  r)e(r)  dr, 


(5.2.20) 


where  H(t)  is  the  Heaviside  step  function.  The  function  ip(t )  represents  the  time  dependent 
behaviour  of  A,  and  therefore  of  y  also,  the  constant  of  proportionality  being  embodied  in  the 
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matrix  D0.  Now  taking,  for  example,  the  Laplace  transform  of  (5.2.20)  which  is  of  the  convolution 
form  we  have 

»(s)  =  D0t?(s)f(s).  (5.2.21) 

Notice  that  any  spatial  variation  of  these  quantities  is  unaffected  by  this  operation. 

Also,  suppose  that  we  know  a  linear  elastic  solution  expressed  in  the  form 


<?e{x)  =  De£e{x) 


(5.2.22) 


which  is  time  independent.  The  correspondence  principle  may  now  be  invoked;  this  consists 
of  replacing  De  in  (5.2.22)  by  DoO(s)  and  treating  <te  and  £e  as  a{s)  and  i(s)  respectively. 
The  result  may  then  be  transformed  inversely  to  yield  the  viscoelastic  solution  corresponding  to 
(5.2.22).  Obviously,  if  the  quantities  in  (5.2.22)  are  time  dependent  this  must  be  incorporated  in 
the  transformation-replacement-inversion  procedure. 

We  consider  the  simple  case  where  the  body  forces  X(t)  and  Y(t)  are  expressed  in  terms  of 
the  single  function 


p(t)  —  ai  +  e  ^‘((^sinwf  +  a2  cos  tot) 
a\  a2to  +  a3(s  +  /?) 

p(!,)  =  7+  (,  +  w+^  ' 


(5.2.23) 

(5.2.24) 


In  the  first  example  we  assume  that  Ar(<)  =  p(t)  and  in  the  second  example  we  assume  that  X(t) 
and  Y(t)  are  given  by  (5.2.11-12)  with  L(t)  =  Lp(t)  and  N(t)  =  Np(t)  where  L  and  N  do  not  vary 
with  time. 

For  simplicity  we  take  the  stress  relaxation  function  < p(t)  of  the  form  (5.1.3)  but  with  only  2 
terms,  i.e. 

<p(t)  =  C0  +  C1e-bt,  C0,  Ci,  6  >  0,  (5.2.25) 

=>  <p(s)  -  go  +  -^"7.  Re(b)  <  Re(s),  (5.2.26) 

s  +  b 

where 

go  '■=  Co  +  Ci  >  0, 
gi  :=  -  bCi  <  0. 


Performing  the  replacement-inversion  procedure  indicated  earlier  we  find  that  for  UE(x,y,t )  and 
vE{x,y,t)  given  by  either  the  quadratic  or  sinusoidal  solution 


u{x,y,t)  _  v(x,y,t) 
ue(x,  y,  t)  vE(x,y,t ) 


ai(feffo  +  gie  ht) 
go(bgo  +  9 1) 


+  ^  -  ur'TZi)  [<J  “ h){a2U + ~  h))‘ 

+  e~ P*  ((/i  —  /?)(6  —  P)  +  w2)  (a2  sin  tot  +  a3  cos  cot) 


- hi 


+  e  —  /t)(a3  sin  cot  —  cos  ut)  L 


(5.2.27) 


where 


bgo  +  g  i 

9o 


(5.2.28) 


Expression  (5.2.27)  furnishes  us  with  the  exact  solution  to  the  corresponding  linear,  non-ageing 
viscoelastic  problem. 
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5.2.3  Viscoelastic  Poisson  effect. 

In  more  general  cases  the  assumption  of  constant  Poisson  ratio  embodied  in  (5.2.18)  will  not  be 
appropriate.  In  this  case  we  must  assume  that  the  time  dependent  behaviour  of  A  and  fi  will  be 
uncorrelated  and  a  deformed  body  may  therefore  exhibit  a  viscoelastic  Poisson  effect,  i.e.  a  time 
dependent  Poisson  ratio. 

Again  using  correspondence  principles,  this  time  under  the  assumption  of  plane  strain  we 
present  an  analytical  solution  to  a  very  simple  problem  which  exhibits  this  time  dependent  Poisson 
effect.  We  use  Cl  as  defined  by  (5.2.1),  with  l  =  1  and  k  =  1/2  and  impose  the  essential  boundary 
condition 

u£(0,  y)  —  0,  u£(l,y)  =  £>  0  (and  constant)  Vy  €  [0, 1] 

»b(z,1/2)=  OVxG  [0,1], 


(5.2.29) 


and  assume  zero  tractions  over  ( x,y )  £  {0  <  x  <  1;  y  =  0,1}  and  identically  zero  body  forces. 
Note  that  the  boundary  condition  (2.2.2)  is  violated  by  this.  In  this  case  the  linear  elastic  solution 
is  given  by 

£11  =  0 

vi  (5.2.30) 


We  may  relate  the  Poisson  ratio  to  the  bulk  modulus,  K,  and  the  shear  modulus,  G,  by  the 
relations  (2.1.5).  Specifically  we  have 


3K  -  2G 

U~  6A+2G 

(5.2.31) 

and  define  the  Lame  coefficients  via 

H-  2G 

(5.2.32) 

To  determine  the  corresponding  viscoelastic  solution  we  assume  the  forms 

K(t)  =  Koexp(-k0t) 

(5.2.33) 

and 

G(t)  =  G0  exp(-got) 

(5.2.34) 

for  the  viscoelastic  bulk  and  shear  moduli,  (the  assumptions  2.1.1  are  satisfied). 

Thinking  again 

of  the  transformation-replacement  procedure  outlined  above,  we  have,  in  the  transform  domain, 

,  3Ko(s  +  </o)  —  2Go(s  +  fco) 

6Ko(s  +  go)  +  2Go(s  +  &o) 

(5.2.35) 

and  thus,  since  l  is 

a  constant: 

with 

A  =  -£(3K0  -  2G0),  B  =  -l(3K0go  -  2G0k0), 

C  =  3/\o  +  4Go,  D  =  ZKogo  +  4Go&o, 

we  may  write 

,_x  ( B/D )  (A/C)  -  (B/D) 

Cl2(s)=  s  +  s  +  (D/C) 

(5.2.36) 

=>  £22(0  =  ^(0 

where 

v(t)  =  No  +  N\  exp (—Nt) 

(5.2.37) 

and 

3/foffo  —  2Goko 

ZKogo  +  AGoko 

s ,  SKogo  —  2Goko  3A'o  —  2Go 

3Aoffo  +  4Gofco  3Ao  +  4Go 
v  3A0  +  4Go 

3Aoffo  +  4Go^o 
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We  may  think  of  i/(t)  as  a  manifestation  of  the  viscoelastic  Poisson  ratio  under  plane  strain,  note 
that  go  =  ko  =>■  £22(0)  =  £22(00)  because  Ni  =  0,  which  is  exactly  the  case  when  it  is  assumed 
that  A  oc  n,  i.e.,  the  Poisson  ratio  is  constant. 

5.3  Numerical  results. 

In  this  section  we  shall  give  graphical  indication  of  the  performance  of  the  algorithms  and 
Pj*fc  by  solving  numerically  the  test  problems  developed  in  the  foregoing  sections.  For  this  we  take, 
in  (5.2.1),  (l,  2k)  =  (1, 1).  Also,  rather  than  calculate  ||u(x,tr)  -  (uh(a:))r||i,n  to  verify  theorems 
4.2.1  and  4.3.1  we  calculate  the  more  interesting  quantity  ||<x(x,<r)  -  (<rfc(x))r||o,n,  which  we  may 
reasonably  expect  to  be  bounded  in  the  same  way  because  the  stress,  <7,  depends  upon  the  space 
derivatives  of  u. 

In  common  with  the  standard  practice  of  computational  finite  element  implementation  we  shall 
perform  all  spatial  integrations  on  a  reference  element  using  a  sufficiently  high  order  quadrature 
rule.  Also  we  take  the  standard  8-noded  quadrilateral  finite  element,  in  which  case  the  basis 
functions  on  this  reference  element  are  of  the  form  Yl!  j  =0  a'j x' ^  >  where  <222  =  0.  The  domain  is 
partioned  into  squares  of  side  length  h  such  that  /i-1  is  an  integer.  Also,  for  simplicity,  we  take 
constant  time  step  k. 

In  order  to  determine  the  order  of  convergence  (i.e.  the  powers  of  h  and  k)  suggested  by  the 
numerical  performance  of  the  algorithms  we  define  ft  through 

R  _  logje2h,2tl  —  log  [eh,fc| 
p:~  log  2 


(5.3.1) 


where  e/,*  is  some  numerically  determined  measure  of  the  error  resulting  from  the  discretization, 
i.e.  /?  is  the  indicated  order  of  convergence  of  the  scheme. 

For  all  of  the  problems  considered  we  assume  the  time  dependent  part  of  the  body  forces 
(5.2.23)  to  take  the  form: 

p(t)  -  1  +  exp(-^f)  sin  .  (5.3.2) 

Also,  for  the  case  of  constant  Poisson  ratio  we  shall  assume  that  v  =  0.4  and  in  (5.2.25) 

ip(t)  —  1  +  exp(— 2f),  (5.3.3) 

whilst  for  the  variable  Poisson  ratio  example  given  in  §5.2.3  we  shall  take  for  (5.2.33,34) 

K(t)  =  8exp(— O.lf),  (5.3.4) 

G(t)  =  4.5  exp(— 0.2t).  (5.3.5) 


We  emphasise  that  the  values  taken  for  the  material  functions  are  purely  arbitrary,  more 
realistic  values  are  considered  in  the  next  section.  Also,  in  order  to  present  the  graphical  results 
in  a  clear  manner  the  solutions  obtained  at  each  time  level  have  been  interpolated  to  produce 
continuous  curves. 

At  the  end  of  §4.3.1  we  remarked  that  the  numerical  performance  of  P*’*  seems  to  indicate 
stability,  in  fact  it  appears  that  if  the  body  forces  and  tractions  tend  to  functions  of  the  space 
variables  alone,  then  for  the  type  of  relaxation  function  given  by  (5.3.3)  the  error  due  to  the  time 
discretization  appears  to  vanish,  leaving  only  the  error  due  to  the  spatial  discretization.  This 
behaviour  does  not  occur  for  P\’k .  Figures  5. 3. 1,3  illustrate  this  by  showing  how  the  particular 
error,  e((l,^),tr)  :=  «i  ((1,  ^),  tr))  -  (uj(l,^))r  varies  through  time.  Also  by  taking  eh,k  to 
be  ||u  —  u^Hoo  (the  maximum  norm  on  TZn)  when  applied  to  the  nodal  solution  vector  for  two 
computations  having  time  steps  k  and  2k  we  show  in  figures  5. 3. 2, 4  the  numerically  suggested 
nodal  convergence  rate,  /?,  for  the  displacements.  For  these  calculations  we  have  used  the  quadratic 
solution  (5.2.7, . . .,  10)  in  which  case  we  expect  the  finite  element  to  capture  the  spatial  behaviour 
exactly. 
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time,  t. 


Figure  5.3.1  Algorithm  P}’0'8,  the  com¬ 
puted  solution  for  (uj (1,  ^))r  and  the  asso¬ 
ciated  error.  The  discretization  error  tends 
to  a  constant  as  the  body  force  tends  to  a 
function  of  x  only. 


time,  t. 

Figure  5.3.2  Algorithm  P}’*,  the  esti¬ 
mated  nodal  convergence  rate  of  the  dis¬ 
placements  with  /?,  given  by  (5.3.1),  ob¬ 
tained  by  considering  k  —  0.8  and  k  =  0.4. 


Figure  5.3.3  Algorithm  Pj'0 '8,  the  com¬ 
puted  solution  for  (uj(l,  |))r  and  the  asso¬ 
ciated  error.  The  discretization  error  tends 
to  zero  as  the  body  force  tends  to  a  func¬ 
tion  of  x  only. 


Figure  5.3.4  Algorithm  P\’k ,  the  esti¬ 
mated  nodal  convergence  rate  of  the  dis¬ 
placements  with  /?,  given  by  (5.3.1),  ob¬ 
tained  by  considering  k  =  0.8  and  k  =  0.4. 


Turning  now  to  the  trigonometric  solution  (5.2.13, ...,  17)  with  £  =  <f>  —  ir  and  L(t),  N(t ) 
given  by  (5.3.2)  above,  figures  5. 3. 5,..., 8  show  the  variation  of  ||<r  —  orh||0in  with  time  and  the 
numerically  suggested  convergence  rate  for  the  two  cases  k  <C  h  and  h  k.  In  each  case  the 
scheme  appears  to  be  0(h 2  +  k2)  as  expected.  The  integral  involved  in  the  calculation  of  ||  •  ||o,n 
was  calculated  numerically  with  an  <D(h 4)  quadrature  rule. 

Finally,  in  figures  5.3.9,10  we  show  the  numerically  predicted  viscoelastic  Poisson  effect.  The 
caption  to  figure  5.3.9  shows  the  values  assumed  for  the  quantities  defined  in  §5.2.3,  and  the  figure 
itself  shows  the  displacement  U2  at  the  point  (1,1),  i.e.  the  top  right  hand  corner  of  the  domain. 
The  accompanying  figure,  5.3.10,  shows  for  interest  only  the  apparent  convergence  rate  of  the 
strain  at  this  point. 

5.4  Comparison  with  a  real  material. 

In  [2]  design  data  for  the  I.C.I.  structural  nylon  66  “Maranyl” (registered  trade  mark),  type  A101 
is  given.  Specifically  we  refer  to  figures  1  and  3  contained  therein.  Figure  3  shows  graphically 
the  response  of  a  test  piece  to  a  constant  imposed  axial  strain  by  plotting  the  isometric  ax¬ 
ial  stress  vs.  time  curves  for  percentage  strains  of  0.5, 1.0, ...,  3.0, 3.5  under  dry  conditions  at 
20° C.  Figure  1  shows  the  tensile  creep  behaviour  resulting  from  constant  imposed  axial  stresses 
of  500,  1000, . . . ,  7000,  7500  p.s.i.  (pounds  per  square  inch).  Thus  figures  1  and  3  represent  collec¬ 
tions  of  axial  creep  and  axial  stress  relaxation  curves  respectively.  Following  the  manufacturers 
recommendation  we  adopt  a  constant  Poisson  ratio  of  0.4. 


time,  t. 


Figure  5.3.5  Algorithm  P}’0'1,  the  com¬ 
puted  value  of  ||<r  — £rh||0  n.  The  discretiza¬ 
tion  error  tends  to  a  constant  as  the  body 
force  tends  to  a  function  of  x  only. 
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time,  t. 

Figure  5.3.6  Algorithm  Pj’4,  the  esti¬ 
mated  convergence  rate  of  ||cr  —  <TA||o,fi  with 
/?,  given  by  (5.3.1),  obtained  by  consider¬ 
ing  h  =  1,  Jfe  =  0.1  and  h  =  0.5 .k  =  0.05. 


•lo- 


time,  t. 

i  i  o 

Figure  5.3.7  Algorithm  P*’  ,  the  com¬ 

puted  value  of  ||<r  —  <Th||0in-  The  discretiza¬ 
tion  error  tends  to  a  constant  as  the  body 
force  tends  to  a  function  of  x  only. 


time,  t. 


Figure  5.3.9  Algorithm  Pj  °’0'8,  the  com¬ 
puted  value  of  £22((^uh(l,  l))r)  and  the 
associated  error,  demonstrating  the  vis¬ 
coelastic  Poisson  effect  discussed  in  §5.2.3, 
in  the  notation  of  that  section  we  have: 
£  =  0.1,  K0  =  8,  Jb0  =  0.1,  Go  =  4.5  and 
go  =  1 


time,  t. 


Figure  5.3.8  Algorithm  Pj,fc,  the  esti¬ 
mated  convergence  rate  of  ||<t  —  <rh||o,n  with 
/?,  given  by  (5.3.1),  obtained  by  consider¬ 
ing  h  —  k  =  1.0  and  h  —  j^,k  =  0.5. 


time,  t. 


Figure  5.3.10  Algorithm  Pj’*,  the  esti¬ 
mated  nodal  convergence  rate  of  £22  —  £22 
at  (1,1)  with  /?,  given  by  (5.3.1),  obtained 
by  considering  k  =  0.8  and  k  =  0.4. 
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In  this  section  we  use  the  0.5%  strain  curve  of  figure  3  to  calculate  an  axial  stress  relaxation 
function  of  the  form  (5.1.3)  by  using  a  least  squares  curve  fitting  routine  (routine  E04GEF  of  the 
NAG  FORTRAN  library,  mkl4).  This  obviously  involves  non-linear  minimisation  and  therefore 
different  minimising  functions  may  result  from  different  starting  points,  we  found  the  following 
coefficients  to  represent  the  stress  relaxation  at  0.5%  strain  adequate: 

C0  =  73  855.4  a0  =  0.0 

Ci  =  155  338.4  a1  =  0.006  139  77  (5.4.1) 

C2  =  173  441.6  a2  =  0.000  181  719 

Where  the  units  of  the  Ci  are  p.s.i,  and  the  a,  are  (hours)-1. 

The  aim  in  this  section  is  to  test  the  predictive  power  of  the  numerical  schemes  by  using  this 
relaxation  function  to  predict  the  creep  curves  given  in  figure  1  of  [2]  (in  fact  for  brevity  we  use 
Pj,fc  only).  The  resulting  behaviour  is  shown  in  figure  5.4.1.  Here  the  data  given  in  [2]  has  been 
reproduced  in  part  by  plotting  the  piecewise  linear  interpolate  to  selected  points  of  figure  1 ,  these 
are  the  solid  lines.  The  calculations  were  performed  under  the  assumption  of  plane  stress. 

The  numerically  computed  results  are  based  on  a  unit  square  domain  discretized  using  a  single 
8-noded  quadrilateral  element  as  in  §5.3  with  the  algorithm  P}’2,  i.e  the  time  step  is  2  hours,  except 
for  the  first  few  steps  where  k  =  0.001  to  capture  the  short  timescale  behaviour.  It  can  be  seen 
that  at  low  stresses  the  algorithm  performs  well  but  the  deterioration  of  accuracy  becomes  very 
marked  as  the  stress  levels  are  increased.  We  conclude  that  the  viscoelastic  relationships  implied 
by  figures  1  and  3  of  [2]  are  more  complicated  than  we  have  allowed  for,  indeed  it  is  most  likely  that 
the  relationship  will  be  non-linear  in  some  way.  The  extension  to  non-linearity  can  occur  in  two 
independent  ways:  firstly  the  assumption  of  infinitesimal  strain  (2.1.3)  may  be  replaced  by  some 
finite  strain  measure,  secondly,  the  linear  constitutive  relationships  (2.1.15,16)  may  be  replaced  by 
some  more  general  form.  Obviously  the  conjunction  of  these  possibilities  gives  the  most  general 
case.  Knauss  and  Emri  [15]  give  an  empirically  determined  formulation  of  non-linear  viscoelasticity 
(also  discussed  in  [16])  based  on  considering  the  free  volume  of  the  material  which  is  used  to  define  a 
reduced  time,  in  this  case  the  infinitesimal  strain  assumption  is  retained.  However,  the  formulation 
of  the  more  general  problem  consistent  with  the  requirements  of  continuum  mechanics  has,  to  our 
knowledge,  not  yet  been  acheived  [17]. 


Figure  5.4.1 


3-11-1992  16:27  P*ge  26 


6  References. 

[1]  Ferry,  J.  D.,  Viscoelastic  Properties  of  Polymers,  John  Wiley  k  Sons,  inc.,  third  edition, 
1980. 

[2] “Maranyl”  nylon,  Technical  Service  Note  N  109,  Moulding  Powders  Group,  I.C.I.  Plastics 

Division,  Welwyn  Garden  City,  Herts,  U.K. 

[3]  Christensen,  R.  M.,  Theory  of  Viscoelasticity,  an  Introduction,  Academic  Press,  London, 
1971. 

[4]  Golden,  J.  M.  &  Graham,  G.  A.  C.,  Boundary  Value  Problems  in  Linear  Viscoelasticity, 
Springer  Verlag,  1988. 

[5]  Johnson,  C.,  Numerical  Solution  of  Partial  Differential  Equations  by  the  Finite  Element 
Method,  Cambridge  University  Press,  1987. 

[6]  Ciarlet,  P.  G.,  The  Finite  Element  Method  for  Elliptic  Problems,  Studies  in  Mathematics 
and  its  applications  4,  North-Holland,  1978. 

[7]  Wheeler,  M.  F.,  A-priori  Hi  Error  Estimates  for  Galerkin  Approximations  to  Partial  Dif¬ 
ferential  Equations,  SIAM  J.  Numer.  Anal.  10  (1973)  pp723-759. 

[8]  Necas,  J.  &  Hlavacek,  I.,  Mathematical  Theory  of  Elastic  and  Elasto-Plastic  bodies, 
an  Introduction,  Studies  in  Applied  Mathematics  3,  Elsevier  Scientific  Publishing  Company, 
Amsterdam,  1981. 

[9]  Fairweather,  G.,  Finite  Element  Galerkin  Methods  for  Differential  Equations,  Lecture  Notes 
in  Pure  and  Applied  Mathematics  vol.  34,  Marcel  Dekker  inc.,  1978. 

[10]  Tricomi,  F.  G.,  Integral  Equations,  Interscience  Publishers,  inc.,  New  York,  fourth  printing 
1967. 

[11]  Whiteman,  J.  R-,  Beagles,  A.  E.  &  Warby,  M.  K.,  Theoretical  and  Practical  Aspects 
of  Finite  Elements  in  the  Context  of  Some  Problems  in  Solid  Mechanincs,  Jber.  d.  Dt. 
Math-Verein,  92  (1990)  pp  77-88. 

[12]  Lees,  M.,  A  Priori  Estimates  for  the  Solutions  of  Difference  Approximations  to  Parabolic 
Partial  Differential  Equations,  Duke  Mathematical  Journal,  27  (1960)  pp  297-311. 

[13]  Linz,  P.,  Analytical  and  Numerical  Methods  for  Volterra  Equations,  Society  for  Industrial 
and  Applied  Mathematics,  1985. 

[14]  Wolkenfelt,  P.  H.  M.,  On  the  Relation  Between  the  Repetition  Factor  and  Numerical 
Stability  of  Direct  Quadrature  Methods  for  Second  Kind  Volterra  Integral  Equations,  SIAM 
J.  Numer.  Anal.  20,  (1983)  pp  1049-1061. 

[15]  Knauss,  W.  G.  &  Emri,  I.,  Volume  Change  and  the  Nonlinearly  Thermo-Viscoelastic 
Constitution  of  Polymers,  Polymer  Engineering  and  Science,  27,  (1987)  pp  86-100. 

[16]  Schapery,  R.  A.,  On  Nonlinear  Viscoelastic  Constitutive  Equations  for  Composite  Materi¬ 
als,  Preprint,  to  appear  in:  Proceedings  of  the  VII  International  Congress  on  Experimental 
Mechanics,  Las  Vegas,  June  8-11,  1992. 

[17]  Whiteman,  J.  R.,  Shaw,  S.  k  Warby,  M.  K.,  Finite  Element  Methods  With  Recovery  for 
Problems  of  Viscoelasticity,  Proceedings  of  International  Conference  on  Numerical  Methods 
in  Engineering  and  Applied  Science,  Concepcion,  Chile,  Nov  1992,  to  appear. 


